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Atomic T Tauri disk winds heated by ambipolar diffusion 
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Abstract. Motivated by recent subarcsecond resolution observations of jets from T Tauri stars, we extend the 
work of ( |Safier 1993^ Jb|) by computing the thermal and ionization structure of self-similar, magnetically-driven, 
atomic disk winds heated by ambipolar diffusion. Improveme nts over his wo rk include: (1) new magnetized cold 
jet solutions consistent with the underlying accretion disk (Ferreira 1997); (2) a more accurate treatment of 



ionization and ion-neutral momentum exchange rates; and (3) predictions for spa tially resolved forbi dden line 
emissi on (maps, lon g-slit spectra, and line ratios), presented in a companion paper, Glarcia et al. (2001). 
As in (Saner 1993c), we obtain jets with a temperature plateau around 10 4 K , but ionization fractions are revised 
downward by a factor of 10-100. This is due to previous omission of thermal speeds in ion-neutral momentum- 
exchange rates and to different jet solutions. The physical origin of the hot temperature plateau is outlined. In 
particular we present three analytical criteria for the presence of a hot plateau, applicable to any given MHD 
wind solution where ambipolar diffusion and adiabatic expansion are the dominant heating and cooling terms. 
We finally show that, for solutions favored by observations, the jet thermal structure remains consistent with 
the usual approximations used for MHD jet calculations (thermalized, perfectly conducting, single hydromagnetic 
cold fluid calculations). 
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1. Introduction 



Norman 1989| ). The main reason lies in the need to pro- 



Progresses in long slit differential astrometry techniques 
and high angular resolution imaging from Adaptive Optics 
and the Hubble Space Telescope have shown that the high 
velocity forbidden emission observed in Classical T Tauri 
Stars (CTTSs) is related to collimated (micro-)jets (eg. 



Solf 1989; |Solf & Bohm 1993; 


Ray et al. 199^; pirth ct al. 


1997]: |Lavalley-Fouquet et al. 2000; 


Dougados et al. 200C; 


Bacciotti ct al. 2000]). Although outflow activity is known 



harbor considerable activity (eg. Mundt fc Eisloffcl 1998 ) 
and present the advantage of not being embedded. It is 
now crtmmonly believed that such jets are magnetically 
self-conhnecl, by a "hoop stress clue to a non-vanisnmg 



duce highly supersonic unidirectional flows. Indeed, this 
requires an acceleration process that is closely related to 
the confining mechanism. The most promising models of 
jet production rely therefore on the presence of large scale 
magnetic fields, extracting energy and mass from a ro- 
tating object. However, we still do not know precisely 
what are the jet driving sources. Moreover, observed jets 
harbor time-dependent features, with time-scales rang- 
ing from tens to thousands of years. Such time-scales are 
much longer than those involving the protostar or the in- 
ner accretion disk. Therefore, although the possibi lity re- 
mains that jets have a non-stationary origin (eg. Ouyecj 



fc Pudritz 1997; Goodson ct al. 199S), only steady-state 



poloidal current (Chan fc Henrikscn 198C; Heyvaerts fc 
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models will be addressed here. 

Stationary stellar wind m odels have been developed 
(eg. Bauty fc Tsinganos 1994 ), however observed correla- 
tions between signatures of accretion and ejection clearly 
show that the disk is an essential ingredient in jet for- 
mation ( Cohen et al. 1989 ; Cabrit et al. 1990| ; Hartigan 
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et al. 1995). Therefore we expect accretion and ejection to 
be interdependent, through the action of magnetic fields. 
There are mainly two classes of stationary magnetized disk 
wind models, depending on the radial extent of the wind- 



heating ( riaficr 1993a| ,|b|) . A further heating scenario (not 
yet explored in the context of MHD jets and only valid 
in some environments) is photoionization from OB stars 
flRcipurth et al. 1998|; [Raga et al. 200$ ft ally fc Reipurth 



producing region in the disk. In the first class (usually 
referred to as "disk winds"), a large scale magnetic field 



threads the disk on a large region ( 


Blandford & Payne 


1982 


3 


'Vardle & Koenigl 1993; 


Ferreira & Pelletier 1993, 


1995 


; 


Li 199E, 1996; Ferreira 1997; 


Krasnopolsky et al. 


1999 


i 


;!asse & Ferreira 20004b 


Vlahakis et al. 2000). Such 



a field is assumed to arise from both advection of interstel- 



lar ma gneti c field and local dynamo generation (Rekowski 
et al. 3000 ). In the second class of models (referred to 
as "X-winds"), only a tiny region around the disk inner 
edge produces a jet (Camenzind 1990 : 5hu et al. 1994 , 



1995| , |1996| ; |Lovelace et al. 1999[) . The magnetic field is 



assumed to originally come from the protostar itself, af- 
ter some eruptive phase that linked the disk inner edge to 
the protostellar magnetosphere. Note that in both mod- 
els, jets extract angular momentum and mass from the 
underlying portion of the disk. However, by construction, 
"disk-winds" are produced from a large spread in radii, 
while "X-winds" arise from a single annulus. Apart from 
distinct disk physics, the difference in size and geometry of 
the ejection regions should also introduce some observable 
jet features. Another scenario has been proposed, where 
the protostar produces a fast collimated jet surrounded 



by a slow uncollimated disk wind or disk corona (Kwan 



& Tademaru 1988, 1995; Kwan 1997), but such a scenario 
still lacks detailed calculations. 

So far, all disc-driven jet calculations used a "cold" 
approximation, ie. negligible thermal pressure gradi- 
ents. Therefore, each magnetic surface is assumed either 
isothermal or adiabatic. But to test which class of models 
is at work in T Tauri stars, reliable observational predic- 
tions must be made and the thermal equilibrium needs 
then to be solved along the flow. Such a difficult task is 
still not addressed in a fully self-consistent way, namely by 
solving together the coupled dynamics and energy equa- 
tions. Thus, no model has been able yet to predict the gas 
excitation needed to generate observational predictions. 

One first possibility is to use a posteriori a simple pa- 
rameterization for the temperature and ionization frac 
tion evolution along the flow 



This was done by Shang 
et al. ( |L998[ ) and |Cabrit et al. (1999| ) for X-winds and disk 



winds respectively. These approaches are able to predict 
the rough jet morphology, but do not provide reliable flux 
and line profile predictions, since the thermal structure 
lacks full physical consistency. 

The second possibility is to solve the thermal evolution 
a posteriori, with the difficulty of identifying the heating 
sources (subject to the constraint of consistency with the 
underlying dynamical solution). Several heating sources 



are i ndeed possible: (1) planar shocks (eg. |Hartigan et al 



1987, 1 



994|); (2) oblique magnetic s hocks in rccollimating 



winds (lOuyed fc Pudritz 1993| , |l994|) ; (3) turbulent mixing 



layers (eg. [Bmcttc ct al. 1999[ ); and (4) current dissipation 
by ion- neutral collisions, referred to as ambipolar diffusion 



2001). Of all these previous mechanisms only ambipolar 
diffusion heating allows "minimal" thermal solutions, in 
the sense that the same physical process — non- vanishing 
currents — is responsible for jet dynamics and heating. As 
a consequence no additional tunable parameter is invoked 



for the thermal description. Furthermore, Bafier (1993b) 
was able to obtain fluxes and profiles in reasonable agree- 
m cnt with obser vations. In this paper, we extend the work 
of Baficr (1993a b|) by (1) using magnetically-driven jet so- 
lutions self-consistently computed with the underlying ac- 
cretion disk, and (2) a more accurate treatment of ioniza- 
tion using the Mappings Ic code and ion-neutral momen- 
tum exchange rates which include the thermal contribu- 
tion. In a companion paper (Garcia et al. 2001, hereafter 



Paper II), we generate predictions for spatially resolved 
orbidden line emission maps, long-slit spectra, and line 
ratios. 

This article is structured as follows: in section 2 we 
introduce the dynamical structure of the disk wind under 
study, and present physical values of the density, veloc- 
ity, magnetic field, and Lorentz force along streamlines ; 
in section 3 we describe the physical processes taken into 
account in the thermal evolution computations, whose re- 
sults are presented and discussed in section 4. Conclusions 
are presented in section 5. Some important derivations, 
dust description and consistency checks of our calculations 
are presented in the appendices. 

2. Dynamical Structure 

2.1. General properties 

A precise disk wind theory must explain how much matter 
is deviated from radial to vertical motion, as well as the 
amount of energy and angular momentum carried away. 
This implies a thorough treatment of both the disk interior 
and its matching with the jets, namely to consider magne- 
tized accretion-ejection structures (hereafter MAES). The 
only way to solve such an entangled problem is to take into 
account all dynamical terms, a task that can be properly 
done within a self-similar framework. 



In this paper, we use the models of Ferreira (1997) 
describing steady-state, axisymmetric MAES under the 
following assumptions: (i) a large scale magnetic field of 
bipolar topology is threading a geometrically thin disk; 
(ii) its ionization is such that MHD applies (neutrals are 
well-coupled to the magnetic field); (iii) some active turbu- 
lence inside the disk produces anomalous diffusion allow- 
ing matter to cross the field lines. Two extra simplifying 
assumptions were used: (iv) jets are assumed to be cold, 
i.e. powered by the magnetic Lorentz force only (the cen- 
trifugal force is due to the Lorentz azimuthal torque), with 
isothermal magnetic surfaces (the midplane temperature 
varying as To oc r _1 ) and (v) jets carry away all disk an- 
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gular momentum. This last assumption has been removed 
only recently by Casse & Ferreira (2000a). 



All solutions obtained so far display the same asymp- 
totic behavior. After an opening of the jet radius leading 
to a very efficient acceleration of the plasma, the jet un- 
dergoes a refocusing towards the axis (recollimation). All 
self-similar solutions arc then terminated, most probably 
producing a shock pomcz de Castro fc Pudritz (199% 



Ouycd & Pudritz (1993). This systematic behavior could 
well be imposed by the self-similar geometry itself and 
not be a general result (Ferreira 1997). Nevertheless, such 
a shock would occur in the asymptotic region, far away 
from the disk. Thus, we can confidently use these solu- 
tions in the acceleration zone, wh ere forbidden emissi on 
lines are believed to be produced ( Hartigan et al. 1995| ). 



2.2. Model parameters 

The isothermal self-similar MAES considered here are de- 
scribed with three free dimensionless local parameters (see 
Ferreira 1997, for more details) and four global quantities: 
(1) the disk aspect ratio 



radius of the MAES is set to Wi = 0.07 AU (typical disk 
corotation radius for a 10 days rotation period): inside 
this region the magnetic field topology could be signif- 
icantly affected by the stellar magnetosphere-disk inter- 
action. The outer radius is kept at zu c = 1AU for con- 
sistency with the one fluid approximation (Appendix |c|) 
and the atomic gas description. Regarding atomic consis- 
tency, |5aficr (1993a ) solved the flow evolution assuming 
initially all H bound in H 2 . He found H 2 to completelly 
dissociate at the wind base, for small m . However, after a 
critical flow line footpoint H 2 would not completelly disso- 
ciate, therefore affecting the thermal solution. This critical 
footpoint was at 3 AU for his MHD solution nearer our 
parameter space. 

We note that our two free parameters are still bounded 
by observational constraints: Mass conservation relates 
the ejection index £ to the accretion/ejection rates ratios, 



2Mj ~ £M acc In — 



(4) 



h(vo) 



-CD- 



The observational estimates for the ratio of mass out- 
flow rate by mass accretion rate are Mj/M acc ~ 0.01 
( Hartigan et al. 1995| ). The uncertainties affecting these 
estimates can be up to a factor of 10 ( pullbring et al 



1998; Lavalley-Fouquct et al. 2000| ). The range of ejec- 



where h(w) is the vertical scale height at the cylindrical 
radius vj; 

(2) the MHD turbulence parameter 



tion indexes considered here (0.005-0.01) is kept compat 
ible with Hartigan's canonical value. The accretion rates 
Mace are also kept free but inside the observed range of 



10~ 8 M a yr ~ 1 to lO _5 M yr _1 in T Tauri stars flHartigan 



— VKh 

where v m is the required turbulent magnetic diffusivity 
and Va the Alfven speed at the disk midplane; this dif- 
fusivity allows matter to cross field lines and therefore to 
accrete towards the central star. It also controls the am- 
plitude of the toroidal field at the disk surface. 
(3) the ejection index 



(2) et al. 1995 



3)- 
0r 



Table |J provides a list of some disk and jet parameters. 
These local parameters were constrained by steady-state 
requirements, namely the smooth crossing of MHD critical 
points. Disk parameters are useful to give us a view of 
the physical conditions inside the disk. Thus, the required 
magnetic field Bo at the disk midplane and at a radial 
distance wo is 



d\nw 



(3) 



S = 0.3 C 



which measures locally the ejection efficiency (£ = in a 
standard accretion disk), but also affects the jet opening 
(a higher £ translates in a lower opening); 

(4) M* the mass of the central protostar; 

(5) vji the inner edge of the MAES; 

(6) w e the outer edge of the MAES, a standard accretion 
disk lying at greater radii. This outer radius is formally 
imposed by the amount of open, large scale magnetic flux 
threading the disk and producing jets; 

(7) M acc , the disk accretion rate fueling the MAES and 
measured at w c . 

For our present study, we keep only £ and M acc as free 
parameters and fix the values of the other five as follows: 
The d isk aspect ratio was measured by Burrows et al. 
(1996D for HH 30 as ~ 0.1 so we fix e = 0.1. The MHD 



VM 



A/a. 



lO- 7 Af yr- 




G 



(5) 



The global energy conservation of a cold MAES writes 



P — IP- 

A arc — l 



jet 



2P r 



ad 



(6) 



where P a cc is the mechanical power liberated by the accre- 
tion flow, P ]ot the total (kinetic, thermal and magnetic) 
power carried away by one jet and P ra d the luminosity ra- 
diated at one surface of the disk. For the solutions used, 
the accretion power is given by 



Pace V 



GM„ AL 



1w\ 



0.177 



turbulence pa rameter is tak en a m = 1 in order to have 
powerful jets ( Ferreira 1997 ). The stellar mass is fixed at 
Af* = 0.5 Mq, typical for T Tauri stars, and the inner 




where Lq is the solar luminosity and the efficiency fac- 
tor rj = {wi/we)^ — {wi/we) depends on both the local 
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Solution 


t 


C 


-fjct / -Prad 




A 


so n 


A 


0.010 


0.729 


1.46 


0.014 


41.6 


50.6 


B 


0.007 


0.690 


1.46 


0.011 


59.4 


52.4 


C 


0.005 


0.627 


1.52 


0.009 


84.2 


55.4 



Table 1. Isothermal MAES parameters. With e — 0.1 and 
a m — 1, the only parameter remaining free is £. Here, the 
magnetic lever arm A, mass load k and initial jet opening 
angle 9q are presented to ease comparison with |5afier 's 
models. However these parameters do not uniquely deter- 
mine the MHD solution. 



Jet density (cm a ) 



Jet velocities (lon/s) 





Ifoignetjc fields (G) 

























* t - 
. t ; 



1 10" 10" 

Lorenta force (dyne/cm*) 




10" 10° 
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Fig. 1. Several wind quantities along a streamline for 
model A (long-dashed line), B (solid) and C (dashed): 
jet nuclei density n, velocity, magnetic field, and Lorentz 
force. For the latter three vectors, poloidal components 
(v p , B p , (J x B) p ) are plotted in black and toroidal com- 
ponents (w^ i?^, (J x B)^) in red. The field line is an- 
chored at vjq = 0.1 AU, around a 1 Mq protostar, with 
an accretion rate M acc = 10~ 6 MQyr _1 . 



ejection efficiency f and the MAES radial extent. Typical 
values for our solutions are r] ~ 0.9. The ratio Pj e t/-frad 
is given in Table |l|. The jet parameters, mass load n 
and magnetic lever arm A , have the same definition as in 
Blandford fc Payne (1982 ). They are gi ven he re to allow 
a comparison with the solutions used in [Bafiei 's work. 



2.3. Physical quantities along streamlines 

In order to obtain a solution for the MAES, a variable sep- 
aration method has been used allowing to transform the 
set of coupled partially differential equations into a set 
of coupled ordinary differential equations (ODEs). Hence, 
the solution in the (zu, z) space is obtained by solving 



for a flow line and then scaling this solution to all space. 
Once a solution is found (for a given set of parameters in 
Section 2.2), the evolution of all wind quantities Q along 
any flow line is given by: 



Q(zu, z) = Qo{w )Q x {x) 



(8) 



where the self-similar variable x — z/vjq measures the 
position along a streamline flowing along a magnetic sur- 
face anchored in the disk at vjq. In particular, the flow 
line shape equation is given by zu(z) — w^{x), where 
the function ^(x) is provided by solving the full dynami- 
cal problem. In Fig. |l| we plot the values of the jet nuclei 
density h and the poloidal and toroidal components of jet 
velocity, magnetic field and Lorentz force. This is done 
for our 3 models, along a streamline with zuq = 0.1 AU, 
M* = 1 Mq and M acc = lO~ 6 M0yr -1 . Values for other 
vjq and M acc can easily be deduced from these plots using 
the following scalings 



n oc 

V oc 

B oc 

J oc 



M acc M* 



Ml C cMl 
MSccMi 



(9) 



The terminal poloidal velocity is v PtOC ~ \J GM^/^-cuq 
so that solutions with smaller opening angles also reach 
smaller terminal velocities, with higher terminal densities. 
The point where reaches a minimum is also the point 
where the jet reaches its maximum width (we call it rec- 
ollimation point), before the jet starts to bend towards 
the axis. The numerical solution becomes unreliable as we 
move away from this point. The MHD solution is stopped 
at the super- Alfvenic point, which is reached nearer for 
higher £. An illustration of the resulting (tu, z) distribu- 
tion of density n and total velocity modulus for model A 
can be found in Fig. 1 of Cabrit et al. (1999). 

3. Flow thermal and ionization processes 

3.1. Main equations 

Under stationarity, the thermal structure of an atomic 
(perfect) gas with density n and temperature T is given 
by the first law of thermodynamics: 



PV ■ v 



V • Uv = T - A, 



(10) 



where P = nkT is the gas pressure, U = ^nkT its internal 
energy density, v the total gas velocity and T and A are 
respectively the heating and cooling rates per unit volume. 
Since during most of the flow the ejected gas expands, we 
call the term A ac jia = -PV • v the adiabatic cooling. 

The gas considered here is composed of electrons, 
ions and neutrals of several atomic species, namely n = 
n c + nl + where the overline stands for a sum over 
all present chemical elements. We then define the density 
of nuclei n = nl + and the electron density n e = / e n. 
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Correspondingly, the total velocity v appearing in Eq. ( |l0| ) 
must be understood as the barycentric velocity. As usual 
in one-fluid approximation, we suppose - and verify it in 
all species well coupled (through collisions), 



C.l 



section 

so that they share the same temperature T. We also as- 
sume that no molecule formation occurs, so that mass con- 
servation requires 



V • hv = 



(11) 



Under stationarity, the gas species ionization state evolves 
according to the rate equations, 



Df A > 
Dt 



v-Vf, 



R 



Element 


Z 


ZdL84 


z d 


H 


1.0 






He 


1.0 (-1) 






C 


3.55(-4) 


3.0 (-4) 


2.17(-4) 


N 


9.33(-5) 




1.39(-5) 


O 


7.41 (-4) 


1.36(-4) 


3.39(-4) 


Ne 


1.23(-4) 






Fe 


3.24(-5) 


3.01(-5) 


3.22(-5) 


Mg 


3.80(-5) 


3.71(-5) 


3.69(-5) 


Si 


3.55(-5) 


3.33(-5) 


3.37(-5) 


S 


1.86(-5) 






Ca 


2.19(-6) 




2.19(-6) 


Ar 


3.63(-6) 




2.43(-6) 


Na 


2.04(-6) 




1.81(-6) 



(12) 
(13) 



Table 2. Abundances by number of various elements with 
respect to Hydrogen. The notation 3.55(-4) means 3.55 x 
1CP 4 . Column Zq gives solar abundances, from Savage 
fe Scmbach (1996| ). Column Zdl84 elemental abundances 



subjected to the elemental conservation constraint, 

Ra> = o. 



(14) 



In the above equations f^t = n^i /h is the population 
fraction of the chemical element A, i times ionized, R^i 
the rate of change of that element and Na is its maxi- 
mum ionization state. Note that R^ is a function of all 
the species densities n^i , the temperature and the radi- 
ation field. In order to obtain the gas temperature and 
ionization state, we must solve energy equation coupled 
to the ionization evolution. This task requires to specify 
the ionization mechanisms as well as the gas heating (r) 
and cooling (A) processes (other than A ac ii a ). 

3.2. Ionization Evolution 
3.2.1. The Mappings Ic code 

We solve the gas ionizat ion state (Eqs. [l2| t o 14 ) using 
the Mappings I c code - Binette et al. (1985 ); Binette fe 
Robinson (1987 ); Fcrruit et al. (1997 ). This code considers 
atomic gas composed by the chemical elements H, He, C, 
N, O, Ne, Fe, Mg, Si, S, Ca, Ar. We also added Na (whose 
ionization evolution is not solved by Mappings Ic), as- 
suming it to be completely ionized in Nan. Hydrogen and 
Helium are treated as five level atoms. 

The rate equations solved by Mappings Ic include 
photoionization, collisional ionization, secondary ioniza- 
tion due to energetic photoelectrons, charge exchange, re- 
combination and dielectronic recombination. This is in 



contrast with Safier , who assumed a fixed ionization frac- 
tion for the heavy elements and solved only for the ion- 
ization evolution of H and He, considering two levels for 
H and only the ground level for He. 

The adopted abund ances are presented in Table 2. 
In contrast with Safier , we take into account heavy el- 
ement depiction onto dust grains (see section 3.2.3 and 
Appendix [bJ) in the dusty region of the wind. 



locked in grains for a MRN-type dust distribution, from 
|Draine fe Lee (1984[ ). Column Zd gives the abundances 
locked in grains in the diffuse clouds toward £ Ophiuchi 
( Savage fe Sembach 1996| )), computed using the solar gas 
phase abundances from the same authors. The depleted 
abundances adopted here arc Z = Zq — Z^. 



3.2.2. Photoionizing radiation Field 

For simplicity, the central source ra diation field is de- 
scribed in exactly the same way as in Safier and we refer 
the reader to the expressions (C1-C10) presented in his 
Appendix C. This radiation field is diluted with distance 
but is also absorbed by intervening wind material ejected 
at smaller radii. 

We treat the radiative transfer as a simple absorption 
of the diluted central source, namely 



47rJ„(r 



L,y(9) 



0) 



4 7rr 2 



(15) 



where J v (r, 9) is the local mean monochromatic intensity 
at a spherical radius r and angle with the disk vertical 9, 
L v (9) is the emitted luminosity of both star and boundary 
layer and t v (r, 9) the optical depth towards the central 
object. 

We now address the question of optical depth. In 
our model, the flow is hollow, starting from a ring lo- 
cated at the inner disk radius Wi and extending to the 
outer radius vj e . The jet inner boundary is therefore ex- 
posed to the central ionizing radiation, which produces 
then a small layer where hydrogen is completely pho- 
toionized. The width Ar of this layer can be computed 
by equating the number of emitted H ionizing photons, 
Q(Ho) = | / L v dvjhv, to the number of recombina- 
tions in this layer, n^aB(T)27rr 2 Ar for our geometry. We 
found that Ar -C r, and thus assume that all photons ca- 
pable of ionizing Hydrogen are exhausted within this thin 
shell. Furthermore, there is presumably matter in the in- 
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ner "hollow" region, so the previous considerations are 
upper limits. 

For the heavy elements, photoionization optical depths 
are negligible, due to the much smaller abundances, and 
are thus ignored. The opacity r„ is assumed to be dom- 
inated by dust absorption (see Appendix ^ for details). 
Dust will influence the ionization structure at the base of 
the flow, where ionization is dominated by heavy elements. 

To summarize, the adopted radiation field is a central 
source absorbed by dust, with a cutoff at and above the 
Hydrogen ionization frequency. 

3.2.3. Dust properties and gas depletion 

Safier showed that if dust exists inside the disk, then 
the wind drag will lift the dust thereby creating a dusty 
wind. Our wind shares the same property. We model the 
dust (Appendix |^) as a mix of graphite and astronomi- 
cal silicate, with a MRN size distribution and use for the 



dust optical properties the tabulated values of Draine & 
Lee (1984| ); |Draine fc Malhotra (1993| ); |Laor fe Draine 
(1993| ) . For simplicity we assumed the dust to be station- 



ary, in thermodynamic equilibrium with the central radi- 
ation field and averaged all dust quantities by the MRN 
size distribution. 

In addition, we take into account depletion of heavy el- 
e ments in to the dust phase. This effect was not considered 
by Safiei . In Table 2 we present the dust phase abundances 
needed to maintain the MRN distribution ( Draine fc Lee ) , 
and our adopted depleted abundances , taken from obser 
rations of diffuse clouds toward £ Ori ( Savage & Scmbach 



1996). These are more realistic, although presenting less 



depletion of carbon than required by MRN. Depletion has 
only a small effect on the calculated wind thermal struc- 
ture, but can be significant when comparing to observed 
line ratios based on depleted elements. 

3.3. Heating & Cooling Mechanisms 
3.3.1. MHD heating 

The dissipation of electric currents J provides a local heat- 
ing term per unit volume Tmhd = J • (E + — X B), where 
E and B are the electric and magnetic fields , v the fluid 
velocity and c the speed of light. In a multi-component 
gas, with electrons and several ion and neutral species, 
the generalized Ohm's law writes 



x B 



2 ±(JxB)xB 



VR J x B 



en c c 



(16) 



where rj is the fluid electrical resistivity, p and p n are 
the total and neutral mass density, rrii n the reduced ion- 
neutral mass, rii the ionic density and Vi n the ion-neutral 
collision frequency. The overline stands for a sum over all 
chemical elements relevant to a given quantity. These last 



quantities depend on the gas ionization state, the temper- 
ature and the momentum exchange rate coefficients. The 
reader is referred to Appendix [A] for the expressions of 
these coefficients and the uncertainties affecting them. 

The first term appearing in the right hand side of the 
generalized Ohm's law is the usual Ohm's term, while the 
second describes the ambipolar diffusion, the third is the 
electric field due to the electron pressure and the last is 
the Hall term. This last effect provides no net dissipation 
in contrary to the other three. It turns out that the dis- 
sipation due to the electronic pressure is quite negligible 
and has been therefore omitted (Appendix |c]) . Thus, the 
MHD dissipation writes 



MHD 



vJ 2 



Ohm 



p£\^\\J x B\ 
P 

r d 



It I , n I ) j Vi a 



(17) 



The first term, Ohmic heating Toiim, arises mainly from 
ion-electron drag. The second term is the ambipolar diffu- 
sion heating Tdrag and is mainly due to ion-neutral drag. 
This last term is the dominant heating mechanism in 
[Safici 's disk wind models, as well as in ours. 

An important difference with Safier is that we take 



into account thermal speeds in ion-neutral momentum ex- 
change rate coefficients. This increases z^ n , and results in 
significantly smaller ionization fractions (Sect. 4.8). 

3.3.2. lonization/recombination cooling 

Both collisional ionization cooling A co i and radiative re- 
combination cooling A roc effects are taken into account by 
Mappings Ic. These terms are given by, 



i<N A 

Acol = ^2 X! R A^A' 
A i=0 
i<N A 

A «c = EE n c n At+ ikTp B (A l 

A i=0 



(18) 



(19) 



where R% is the collisional ionization rate, 1^ the ion- 
ization energy and /3b (^4*) is the case B radiative recom- 
bination rate to the ionization state A 1 . 

These ionizat ioiV rec ombination effects, taken into ac- 
count in part by Safier, are in general smaller than adia- 
batic and line cooling. 



3.3.3. Photoionization heating and radiative cooling 

Photoionization by the radiation field, not taken into ac- 



count by Safier, provides an extra source of heating Tp. 



This term, which is also computed by Mappings Ic, is 
given by 



^^af(hv-hvf)dv (20) 
hv 



i<N A 
A i=0 

where v§ is the threshold frequency for ionization of the 
chemical element A, i times ionized, Air J„ is the radiation 
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flux of the central source (described in section 3.2.2) and 



the photoionization cross-section. We found it to be 
the dominant heating source at the base of the flow, at 
the inner radii and for high accretion rates. 

Collisionally excited lines provide a very efficient way 
to cool the gas, thanks to an extensive set of resonance and 
inter-combination lines, as well as forbidden lines. This ra- 
diative cooling A rac i is computed by Mappings Ic by solv- 
ing for each atom the local statistical equilibrium, and 
will allow us to compute emission maps and line profiles 



for comparison with observations (see Paper II). We in- 
clude cooling by hydrogen lines, A ra d(H), in particular Ha, 
which could not be computed by Saficr| (two-level atom 
description) . 

3.3.4. Other minor heating/cooling mechanisms 

Several processes, also computed by Mappings Ic, ap- 
peared to be very small and not affecting the jet ther- 
mal structure. We just cite them here for completeness: 
free-free cooling and heating, two-photon continuum and 
Compton scattering. 

We ignored thermal conduction, which could be rele- 
vant along flow lines, the magnetic field reducing the gas 
thermal conductibility in any other direction. Also ignored 
was gas cooling by dust grains and heating by cosmic rays. 
We checked a posteriori that these three terms indeed have 



a negligible contribution (see Appendix C.2) 



3.4. Numerical resolution 

In our study, flow thermodynamics are decoupled from the 
dynamics - cold jet approximation. The previous equa- 
tions (|l0|) to ( |l4| ) can then be integrated for a given flow 
pattern. The dynamical quantities (density, velocity and 
magnetic fields) are given by the cold MHD solutions pre- 
sented in Section |^. For the steady-state, axisymmetric, 
self-similar MHD winds under study, any total derivative 
writes 

D , „ v, d , 

— = (t> . V) = — — (21) 

where v z is the vertical velocity and \ — z/zuq is the self- 
similar variable that measures the position along a flow 
line anchored at wq. With this in hand, and the mass 
conservation condition for an atomic wind (Eq. pi] ), the 
energy equation Eq. [lO] becomes an ODE along the flow 
line: 

dlnT _ 2d\nn d + ^ ^ . T-A 
dx 3 d\ d\ 



1 ro 



2 dhin 

3 dx 



1 - 



r 



MHD 



+ (T - A) Map 



A 



adia 



(22) 



The term d\n(l + f e )/dx (due to variations in internal en- 
ergy density) is in fact negligible and has not been imple- 
mented for computational simplicity (see Appendix |C.2|) . 
The term (r - A) Map = T P - A rad - A col - A rcc is pro- 
vided by Mappings Ic. The wind thermal structure is 



computed by integrating along a flow line the energy equa- 
tion ( p2| ) coupled to the set of electronic population equa- 
tions (Eqs. (Ill) to (|TT 



3.4.1. Initial Conditions 



The integration of the set equations (|12|) to ( |14| ) and (22) 
along the flow is an initial value problem. Thus, some way 
to estimate the initial temperature and populations must 
be devised. All calculations start at the slow-magnetosonic 
(SM) point, which is roughly at two scale heights above 
the disk midplane (for the solutions used here). 

To estimate the initial temperature, Safier equated the 
poloidal flow speed at the SM point to the sound speed. 
Although this estimation agrees with cold flow theory, it 
is inconsistent with the energy equation which is used fur- 
ther up in the jet. Our approach was then to compute the 
initial temperature and ionic populations assuming that 



DT 
~Dt 



= and 



Dn A i 



Dt 



= 



(23) 



is fulfilled at x — Xs- We thus assume for convenience that, 
at the base of the jet, there is no strong variations neither 
in temperature nor in ionization fractions. Physically this 
means that the gas is in ionization equilibrium, at the 
obtained temperature, with the incoming radiation field. 
The temperature thus obtained is always smaller than that 
provided by the SM speed. This is due to the large opening 
of the magnetic surfaces, providing a dominant adiabatic 
cooling over all heating processes. For numerical reasons 
the minimum possible initial temperature was set to 50 K. 
It is noteworthy that the exact value of the initial tem- 
perature only affects the base o f th e flow, below a few 
thousand degrees (see Appendix C.3). These regions are 
too cold to contribute significantly to optical line emission, 
leaving observational predictions unaffected. 

The initial populations are computed by Mappings Ic 
assuming ionization equilibrium with the incoming radi- 
ation field. However for high accretion rates and for the 
outer zones of the wind, dust opacity and inclination ef- 
fects shield completely the ionizing radiation field. The 
temperature is too low for collisional ionization to be ef- 
fective. The ionization fraction thus reaches our prescribed 
minimum - all Na is in the form Nail (Table 2) and all 
the other elements (computed by Mappings Ic) neutral. 
However, soon the gas flow gains height and the ioniza- 
tion field is strong enough such that the ionization is self- 
consistently computed by Mappings Ic. 

3.4.2. Integration procedure 

After obtaining an initial temperature and ionization state 
for the gas we proceed by integrating the system of equa- 
tions. In practice the ionization evolution is computed by 
Mappings Ic and separately we integrate Eq. (^2j) with 
a Runge-Kutta type algorithm (Press et al., 1988). We 
maintain both the populations and Mappings Ic cool- 
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ing/hcating rates per n 2 fixed during each temperature in- 
tegrating spatial step. After we call Mappings Ic to evolve 
the populations and rates, at the new temperature, during 
the time taken by the fluid to move the spatial step. This 
step is such, that the RK integration has a numerical ac- 
curacy of 10 -6 and, that the newly computed temperature 
varies by less than a factor of 10 -4 . Such a small varia- 
tion in temperature allows us to assume constant rates 
and populations while solving the energy equation. We 
checked a few integrations by redoing them at half the 
step used and found that the error in the ionic fraction 
is < 10~ 3 in the jet, and < 10~ 2 in the recollimation 
zone; the temperature precision being roughly a few times 
better. This ensures an intrinsic numerical precision com- 
fortably below the accuracy of the atomic data and the 
(cr H H + v) collision cross-sections which, coupled to abun- 
dance incertitudes, are the main limitating factors. Details 
on the actual numerical procedure used by Mappings Ic 
to compute the non-equilibrium gas evolution are given in 
Binctte et al. (1985[). 



4. Thermal structure results 

In this section we present the calculated thermal and ion- 
ization structure along wind flow lines, discuss the physical 
origin of the temperature plateau and its connection with 
the underlying MHD solution, discuss the effect of vari- 
ous key model parameters and finally compare our results 
with those found by Safier. The parameters spanned for 
the calculation of the thermal solutions are the wind ejec- 
tion index £ describing the flow line geometry, the mass 
accretion rate M acc and the cylindrical radius tuq where 
the field is anchored in the disk. 



4.1. Temperature evolution 

In Figure |5[ solid curves present the out of ionization equi- 
librium evolution of temperature, electronic density, and 
proton fraction along flow lines with wq = 0.1 and 1 AU, 
as a function of x = z/tuq, for accretion rates ranging from 
10 -8 to 10 _5 Mq yr _1 . For comparison purposes, dashed 
curves plot the same quantities calculated assuming ion- 
ization equilibrium at the local temperature and radiation 
field. For compactness we present only these detailed re- 
sults for our model B, with an intermediate ejection index 
£ = 0.007. We divide the flow in three regions: the base, 
the jet and the recollimation zone. These regions are sep- 
arated by the Alfvcn point and the recollimation point 
(where the axial distance reaches its maximum). We only 
present the initial part of the recollimation zone here, be- 
cause the dynamical solution is less reliable further out, 
where gas pressure is increased by compression and may 
not be negligible anymore. Note that the recollimation 
zone was not yet reache d over the scales of interest in the 
solutions used by Safier . 

The gas temperature increases steeply at the wind base 
(after an initial cooling phase for high M acc > 1O -6 M 




1 10" 10" 

X ::::: z/^o 

Fig. 2. Several wind quantities versus x — z/ujq for 
model B. The out of ionization equilibrium calculations 
are the solid curves and, for comparison, the ionization 
equilibrium are the dashed ones. The vertical dotted lines 
mark the Alfven point and recollimation point. Top: 
Temperature, Middle: electronic density n c , Bottom: 
proton fraction / p = n(H + )/nn- The accretion rate M acc 
increases in the direction of the arrow from 10~ 8 to 
1O~ 5 M yr _1 in factors of 10. 



around T ~ 1 — 3 x 1Q 4 K, before increasing again af- 
ter the recollimation point through compressive heating. 
The plateau is reached further out for larger accretion 
rates and larger vuq. Its temperature decreases with in- 
creasing M acc . The temperature plateau and its behavior 
with Macc were first identified by Safiei in his wind solu- 
tions. We will discuss in Section 4.4 why they represent a 
robust property of magnetically-driven disk winds heated 
by ambipolar diffusion. 



yr 1 ). It then stabilizes in a hot temperature plateau 



4.2. Ionization and electronic density 

The bottom panels of Fig. || plot the proton fraction 
f p = n(R + )/nn along the flow lines. It rises steeply with 
wind temperature through collisional ionization, reach- 
ing a value ~ 10~ 4 at the beginning of the tempera- 
ture plateau. Beyond this point, it continues to increase 
but starts to "lag behind" the ionization equilibrium cal- 
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culations (dashed curves): the density decline in the ex- 
panding wind increases the ionization and recombination 
timescales. Eventually, for \ ^ 100, density is so low 
that these timescales become longer than the dynamical 
ones, and the proton fraction becomes completely "frozen- 
in" at a constant level, typically a factor 2-3 below the 
value found in ionization equilibrium calculations (dashed 
curves) . 

The electron density (n e ) evolution is shown in the 
middle panels of Fig. ^| In the jet region, where f p is 
roughly constant, the dominant decreasing pattern with \ 
is set by the wind density evolution as the gas speeds up 
and expands. Similarly, the rise in n e in the recollimation 
zone is due to gas compression. A remarkable result is 
that, as long as ionization is dominated by hydrogen (i.e. 



ss a = o.i AV 



a B = 1 AU 

1 " 1 ,M H 



> 10" 



is not highly dependent of M acc , increasing 



by a factor of 10 only over three orders in magnitude in 
accretion rate. This indicates a roughly inverse scaling of 
f p with Macc (b ottom panels of Fig. ||), a property already 



found by Safier which we will discuss in more detail later. 

In regions at the wind base where f p < 10~ 4 , variations 
of n e are linked to the detailed photoionization of heavy el- 
ements which are then the dominant electron donors. The 
respective contributions of various ionized heavy atoms 
to the electronic fraction / e is illustrated in Fig. || for 
M acc = KT 6 M Q yr _1 . While On and Nil are strongly 
coupled to hydrogen collisional ionization through charge 
exchange reactions, the other elements are dominated by 
photoionization. The sharp discontinuity in C II and Na II 
at the wind base for zuq = 0.1 AU is caused by the cross- 
ing of the dust sublimation surface by the streamline (see 
Appendix |b|) . Inside the surface we are in the dust sub- 
limation zone where heavy atoms are consequently not 
depleted onto grains and hence have a higher abundance. 
In contrast, for tzjq = 1 AU, the flow starts already outside 
the sublimation radius, in a region well-shielded from the 
UV flux of the boundary-layer, where only Na is ionized. 
Extinction progressively decreases as material is lifted 
above the disk plane and sulfur, then carbon, also become 
completely photoionized. 



4.3. Heating and cooling processes 

The heating and cooling terms along the streamlines for 
our out of equilibrium calculations are plotted in Fig. 4.3 
for wq — 0.1 and 1 AU, and for two values of M acc = 10~ y 
and 10~ 7 Mq yr" 1 . 

Before the recollimation point, the main cooling pro- 
cess throughout the flow is adiabatic cooling A a di a , al- 
though Hydrogen line cooling A ra( j(H) is definitely not 
negligible. The main heating process is ambipolar diffusion 
Tdrag- The only exception occurs at the wind base for small 




Fig. 3. Ion abundances with respect to hydrogen (f^ = 
n^i jh and thus depends also on the abundances) along 
the flow line versus x = z/^o for model B in out of ion- 
ization equilibrium, with M acc = 10 -6 MQyr -1 . The jump 
at x ~ 0.5 is due to depletion as the gas enters the subli- 
mation surface. 



decays very fast due to the combined effects of radiation 
dilution, dust opacity, depletion of heavy atoms in the dust 
phase, and the decrease in gas density. At the same time, 
the latter two effects make Tdrag rise and become quickly 
the dominant heating term. In the recollimation zone, the 
main cooling process is Hydrogen line cooling A ra d(H), 
and the main heating term is compression heating (A a dia 
is negative). 



A striking result in Fig. 4.3, also found by Bafier, is that 
a close match is quickly established along each streamline 
between A a d; a and r<j ra g! and is maintained until the rec- 
ollimation region. The value of x where this balance is 
established is also where the temperature plateau starts. 
We will demonstrate below why this is so for the class of 
MHD wind solutions considered here. 



4.4. Physical origin of the temperature plateau 

The existence of a hot temperature plateau where A a di a 
exactly balances Tdrag is the most remarkable and robust 
property of magnetically-driven disk winds heated by am- 
bipolar diffusion. Furthermore, it occurs throughout sev- 
eral decades along the flow including the zone of the jet 
that current observations are able to spatially resolve. 

In this section, we explore in detail which generic prop- 
erties of our MHD solution allow a temperature plateau 
at ~ 10 4 K to be reached, and why this equilibrium may 
not be reached for other MHD wind solutions. 



vjq < 0.1 AU and large M acc > 10~ b M Q yr , where pho- 4.4.1. Context 
toionization heating Tp initially dominates. Under such 
conditions, ambipolar diffusion heating is low due to the 
high ion density, which couples them to neutrals and re- 
duces the drift responsible for drag heating. However, Tp 



First, we note that the energy equation (Eq. [2 4 ) in the 
region where drag heating and adiabatic cooling are the 
dominant terms (which includes the plateau region) can 
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Fig. 4. Heating and cooling processes (in erg s _1 cm~ 3 ) along the flow line versus x = z/vjq for model B. Top figures 
for M acc = lO^Moyr" 1 and bottom ones for M acc = 10 _7 Moyr _1 . Ambipolar heating and adiabatic cooling appear 
to be the dominant terms, although Hydrogen line cooling cannot be neglected for the inner streamlines. 



be cast in the simplified form: 



d\nT 
where: 

^(xr 1 



2 din? 



rdra 

a" 



aclia. 



3 dlnx 

\F(T) /' 



- 1 



V 3 



2 dlnj 



3 d\nx 



(24) 



(25) 



remains positive before recollimation and depends only on 
the MHD solution, 



G( X ) = - 



J x B\ 



Mi 



n 2 (v ■ V)n 



M acc w 



(26) 



is a positive function, before recollimation, that depends 
only on the MHD wind solution, and 



F(T) = kT(l + f e )m in fifn (<r in v) x (JL) 



(27) 



also positive, depends only on the local temperature and 
ionization state of the gas. The functions G and F separate 
the contributions of the MHD dynamics and ionization 
processes in the final thermal solution. 



The function 5 is roughly constant and around unity 
before recollimation, it diverges at the recollimation point 
and becomes negative after it. Throughout the plateau 
■) ■■ 1 . 

The "wind function" G is plotted in the center panel of 
Fig. ^|for our 3 solutions. It rises by 5 orders of magnitude 
at the wind base and then stabilizes in the jet region (until 
it diverges to infinity near the recollimation point). The 
physical reason for its behavior is better seen if we note 
that the main force driving the flow is the Lorentz force: 



G 



Dv 



Dt 



\Dt) 



(28) 



For an expanding and accelerating flow (Dp/ Dt) -1 is an 
increasing function. At the wind base the Lorentz force 
accelerates the gas thus causing a fast increase in G. Once 
the Alfven point is reached, the acceleration is smaller and 
Dv/Dt decreases. However this decrease is not so abrupt 
as in the case of a spherical wind, because the Lorentz 
force is still at work, both accelerating and collimating the 
flow. This collimation in turn reduces the rate of increase 
of (Dp/Dt)^ 1 . The stabilization of G observed after the 
Alfven point is thus closely linked to the jet dynamics. 

The "ionization function" F is in general a rising 
function of T and is plotted in the left panel of Fig. || 
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Fig. 5. Left: Function F(T) in erg g cm 3 s _1 versus temperature assuming local ionization equilibrium and an 
ionization flux that ionizes only all Na and all C. Center: Function G(x) in erg g cm 3 s _1 for models A, B, C (bottom 
to top), M acc = lO _6 M0yr _1 and vjq — 0.1AU. Right: Temperature for model B from the complete calculations in 
ionization equilibrium (dashed), and assuming T — Tq as given by eq. [30] ( dash-dotted), wq = 0.1 AU and accretion 
rates M acc are 1CT 8 to 10 -5 M Q yr~ 4 , from top to bottom. 



under the approximation of local ionization equilibrium. 
Two regimes are present: In the low temperature regime, 
/i ^ fp is dominated by the abundance of photoionized 
heavy elements and F(T) cx T/j increases linearly with 
T, for fixed /j. The effect of the UV flux in this region 
is to shift vertically F(T): for a low UV flux regime only 
Na is ionized and /j ~ /(Nan); for a high UV flux regime 
were Carbon is fully ionized, /j ~ /(C n). In the high tem- 
perature regime (T > 8000 K) where hydrogen collisional 
ionization dominates, fi ~ f p , and F(T) cx T f p becomes 
a steeply rising function of temperature, until hydrogen is 
fully ionized around T ~ 2 x 10 4 K. The following second 
rise in F(T) is due to Helium collisional ionization. As we 
go out of perfect local ionization equilibrium the effect is 
to decrease the slope of F(T) in the region where H ion- 
ization dominates. In the extreme situation of ionization 
freezing, F(T) becomes linear again as in the photoionized 
region: F{T) cx T/ Pifrcczcd . 

4.4.2. Conditions needed for a hot plateau 

The plateau is simply a region where the temperature does 
not vary much, 



d\nT 



with 



N « l. 



dlnx 

Naively, temperatures Tq(x) defined by, 
F(T e ) = G( X ) «=► r drag = A adia , 



(29) 



(30) 



will zero the right hand side of Eq. |2j and thus satisfy the 
plateau condition. This equality is the first constraint on 
the wind functiond G, because there must exist a temper- 
ature T e such that the equality holds. However this condi- 
tion is not sufficient. Indeed the above equality describes 



a curve[| T e (x) which must be flat in order to satisfy the 
plateau condition (Eq. ^). Therefore the requirement of 
a flat Tq translates in a second constraint on the varia- 
tion of G with respect to F. This constraint is obtained 
by differentiating Eq. pO. We obtain 



d\uG d\uF 



d In x d In T 



x e 



and after using (Eq. E3): 



rflnG 
dlnx 



< 



,d\nF 



d\nT 



(x)l 



(31) 



(32) 



T=T«= 



Thus only winds where the wind function G varies much 
slower than the ionization function F will produce a 
plateau. This is fulfilled for our models: Below the Alfvcn 
surface, G varies a lot, but collisional H ionization is suffi- 
ciently close to ionization equilibrium that F(T) still rises 
steeply around 10 4 K (Fig. |J). For our numerical values of 
G, within our range of M acc and wo, we have T e ~ 10 4 K 
and thus |dlnG/dlnx| < \dlnF/dhxT\. Further out, 
where ionization is frozen out, we have dhiF/dhiT = 1 
(because F cx T/ p f reezed ) but it turns out that in this re- 
gion G is a slowly varying function of Xi an d thus we still 
have IdlnG/dlnxl < \d\n F / d\nT\. 

Finally, a third constraint is that the flow must quickly 
reach the plateau solution T{x) = Tq(x) and tend to 
maintain this equilibrium. Let us assume that T = Tq 
is fulfilled at x = xo, what will be the temperature at 
X — Xo(l + ? Letting T — Tq(1 + $) and assuming 
eCl. Eq. gives us d — exp(— ax), which provides an 
exponential convergence towards T = Tq as long as 



1 dlnF 



a = — - 



5 rflnT 



> 



(33) 



T=T fi 



1 This is only true if F is a monotonous function of T. As 
can be seen from Figure H this is true for almost all its domain. 
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Note that a depends mainly on the MHD solution (inde- 
pendent of M^cc/wq) and that the steeper the function 
F, the faster the convergence. The above criterion is al- 
ways fulfilled in the expanding region of our atomic wind 
solutions, where 5 > and F increases with temperature. 
The physical reason for the convergence can be easily un- 
derstood in the following way: it can be readily seen that 
if at a given point T > T e (x), then G(x)/F(T) < 1 and 
the gas will cool (cf. Eq. [24| with S > 0). Conversely, if 
T < Tq(x), the gas will heat up. Thus, for S > 0, the fact 
that F(T) is a rising function introduces a feedback that 
brings and maintains the temperature close to its local 
equilibrium value Tq{x), and A a di a close to Fdrag- 

We conclude that three analytical criteria must be met 
by any MHD wind dominated by ambipolar diffusion heat- 
ing and adiabatic cooling, in order to converge to a hot 
temperature plateau: 

(1) Equilibrium: the wind function G must be such that 
F(T) = G(x) is possible around T ~ 10 4 K; 

(2) Small temperature variation: the wind function G{x) 
must vary slower than the ionization function F(T) such 
that \d\nG/d\nx\ < \d\nF/d\nT\: (i) the wind must be 
in ionization equilibrium, or near it, in regions where G is 
a fast function of Xi (ii) once we have ionization freezing, 
G must vary slowly, with \dhiG/dhix\ <C 1; 

(3) Convergence: a > 0, i.e. | (d In n/d In x) x 
(dlnF/dlnT)|r=T e < 0, which is always verified for an 
atomic and expanding wind. 



4.4.3. Comments on other MHD winds 

Not all types of MHD wind solutions will verify our first 
criterion. Physically, the large values of G{x) observed in 
our solutions indicates that there is still a non-negligible 
Lorentz force after the Alfven surface. In this region 
(which we call the jet) the Lorentz force is dominated 
by its poloidal component which both collimates and ac- 
celerates the gas (Fig. p. The gas acceleration translates 
in a further decrease in density contributing to a further 
increase in G. Models that provide most of the flow accel- 
eration before the Alfven surface might turn out to have 
a lower wind function G(x), not numerically compatible 
with the steep portion of the ionization function F(T). 
These models would not establish a temperature plateau 
around 10 4 K by ambipolar diffusion heating. They would 
either stabilize on a lower temperature plateau (on the lin- 
ear part of F(T)) if our second criterion is verified, or con- 
tinue to cool if G varies too fast for the second criterion to 
hold. This is the case in particular for the analytical wind 



models considered by Ruden et al. (199C), where the drag 



force was computed a posteriori from a prescribed velocity 
field. The G function for their parameter space (Table 3 
of Ruden et al.) peaks at ~ 10 -48 erg g cm 3 s _1 at ~ 3i?* 
and then rapidly decreases as G oc r^ 1 for higher radii. 
This translates into a cooling wind without a plateau. 



(M at 



(M„, 



Fig. 6. Verification of the plateau scalings (Eq. pj). Left: 
we plot the measured T e / P ,e versus (M acc ra7o) for all 
models in out of ionization equilibrium. The evolution is 
linear except for the very edges of the (Ma^roo)" 1 do- 
main. This is due to the failure of our assumptions: in 
the lower edge f c < 10~ 4 and thus it isn't dominated by 
H; in the upper edge / c > 0.1 and thus the small ion- 
ization fraction approximation fails. The crosses are from 
[Safici ; because of a smaller momentum transfer rate coef- 
ficient they are quite above model C. Right: we plot the 
measured / PiQ versus (M acc rc7o) _1 for Model B in out of 
ionization equilibrium (solid) and ionization equilibrium 
(dashed). Note that all the variation is absorbed by / Pi q 
and thus / p ,e oc (M acc tuo)" 1 . A straight line is also plot 
for comparison. 



4.5. Scalings of plateau parameters with M acc and wq 

The balance between drag heating and adiabatic cooling 
(Eq. ^o|) can further be used to understand the scalings 
of the plateau temperature Tq and proton fraction / Pj q 
with the accretion rate M acc and flow line footpoint radius 
(wq). In the plateau region, ionization is intermediate, i.e., 
sufficiently high to be dominated by protons but with most 
of the Hydrogen neutral. Under these conditions we have 
F(T) oc Tq/ P) q. On the other hand, self-similar disk wind 
models display G(x) oc (A/ acc tuo) _1 (Eq. ||^). Therefore, 
we expect 



T e f P 



1 



M acc ru a 



(34) 



This behavior is verified in the left panel of Fig. [j] for our 
out-of-equilibrium results for the 3 wind solutions. 

To predict how much of this scaling will be absorbed by 
Tq and how much by / p . e , [Bafici considered the ionization 
equilibrium approximation: For the temperature range of 
the plateau (T ~ 10 4 K) / p ~ /; is a very fast varying 
function of T that can be approximated as / p oc T a with 
a > 1. One predicts that T e oc (M acc ra )~ 1/(a+1) while 
/p, e oc (Af acc n7 )- a /( a+1 > ~ (M acc wo)" 1 . Hence, the in- 
verse scaling with (M acc n7 ) should be mostly absorbed 
by / Pl e, while the plateau temperature is only weakly de- 
pendent on these parameters. This is verified in the right 
panel of Fig. g, where f Py Q in ionization equilibrium is 
plotted as a function of (M acc too)™ 1 . The predicted scal- 
ing is indeed closely followed. 
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Let us now turn back to the actual out-of-cquilibrium 
calculations. At the base of the flow we find that the wind 
evolves roughly in ionization equilibrium (see Figure ^), 
however at a certain point the ionization fraction freezes 
at values that are near those of the ionization equilib- 
rium zone. This effect implies that the ionization fraction 
should roughly scale as the ionization equilibrium values 
at the upper wind base. This in indeed observed in Fig. ^. 
This memory of the ionization equilibrium values by f p 
(as observed in the solar wind by Owocki et al. (1983)) 
is the reason why the scalings of / P) q with (M acc -n7o) re- 
main correct. We computed for our solutions the scalings 
and found for model B: / p , e oc M~°- 76 , / p>e oc w " - 83 , 
Tq oc M~ c ® 13 and no dependence of T on wq, confirming 
the memory effect on the ionization fraction only. 

Finally, we note that for accretion rates in excess of 
a few times 10 -5 M© yr _1 , the hot plateau should not 
be present anymore: Because of its inverse scaling with 



the wind function G remains below 10 



-47 



erg g 



ear low-temperature part of F(T) where /; ~ /(Cll) (see 
Figure ||). These colder jets will presumably be partly 
molecular. Interestingly, molecular jets have only been ob- 
served so far in embedded protostars with high accretion 
rates (e.g. Gueth & Guilloteau 1999). 

4.6. Effect of the ejection index £ 

The importance of the underlying MHD solution is illus- 
trated in figure |[ The ejection index £ is directly linked 
to the mass loaded in the jet (k, see Tab. |l| and Ferreira 
1997). Thus a higher £ translates in an stronger adiabatic 
cooling because more mass is being ejected. The ambipolar 
diffusion heating is less sensitive to the ejection index, be- 
cause the density increase is balanced by a stronger mag- 
netic field. Hence, the wind function G(\) decreases with 
increasing £. As a result, the plateau temperature and ion- 
ization fraction also decrease (see Eq. |34| ). 

In figure [7] we summarize our results for the three mod- 
els by plotting the plateau f Pi Q versus Tq, for several M acc 
(values of wq are connected together). In this plane, our 
MHD solutions lie in a well-defined "strip" located below 
the ionization equilibrium curve, between the two dotted 
curves. For a given model, as M acc increases, the plateau 
ionization fraction and the temperature both decrease, as 
expected from the scalings discussed above, moving the 
model to the lower-left of the strip. Increasing the ejec- 
tion index decreases G{x) 1 and it can be seen that this 
has a similar effect as increasing M acc (Eq. |3|). 

4. 7. Depletion effects on the thermal structure 

In our calculations we take into account depletion of heavy 
species into the dust phase. We ran our model with and 
without depletion and found these effects to be minor. 
Changes are only found when f p < 10~ 4 . The temperature 
without depletion is slightly reduced (the higher ioniza- 




T (K) 

Fig. 7. Out of ionization equilibrium evolution of wind 
quantities in the plateau. Points are for all flow lines (wq = 
0.07 x 1.3 1 AU and i — 0,1,. ..,10) and accretion rates 
(Af acc = KT* Moyr- 1 with i = 5,6,7,8). We plot, / p>e 
versus Tq for all models C (red), B (black), A (blue) and in 
red the values expected from ionization equilibrium. The 
dashed/dotted lines are for models without depletion, the 
solid lines are for models with depletion. The accretion 
rates increase in the direction top right to bottom left. 
The thick solid curve traces / p in ionization equilibrium, 
while the two dotted lines embrace the locus of our MHD 
solutions. 



tion fraction reduces Tdrag) and as a consequence / p is also 
smaller. Normally these changes affect only the wind base, 
as the temperature increases f p dominates the ioniza- 
tion and we obtain the same results for the plateau zone. 
However for high accretion rates (M acc = 10 -5 MQyr -1 ) in 
the outer wind zone (large wq) we still have f p < 10 -4 for 
the plateau and thus the temperature without depletion 
is reduced there. 



4.8. Differences with the results of Safiei 



T he mos t striking difference between our results and those 
of Saficr is an ionization fraction 10 to 100 times smaller. 
This difference is mainly due to both different (ctjj h+ w ) 
momentum transfer rate coefficients and dynamical MHD 
wind models. 

The critical importance of the momentum transfer rate 
coefficient ( (ct h h + v) ) for the plateau ionization fraction 
can be seen by repeating the reasoning in the previous sec- 
tion but including the momentum transfer rate coefficient 
in the scalings. We thus obtain 



T e fp,e « 



1 



(<T H R +v)M &cc W 



(35) 



This shows that, because the freezing of the ionization 
fraction is correlated to the ionization fraction at the base 
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Fig. 8. (<t h H +v) in cm 3 s -1 using Draine expres- 
sion (soli d), using Geiss & Buergi expression (see 
Appendix |A.2| ) (dash) and Draine expression but ignor- 
ing the thermal velocity (dot). 



of the wind (which is in ionization equilibrium), f p will 
scale with the momentum transfer rate coefficient value. 
This means that if the momentum transfer rate coeffi- 
cient is larger, there is a better coupling between ions and 
neutrals and hence a sma ller dr ag heating. For the cal- 
culation of the (<7h h+ v ) i Safiei ignored the contribution 
of the thermal velocity in the collisional relative velocity. 
This considerably reduces (cjh n+ v ) an d thus, increases 
/ p . In figure |§| we plot the corresponding momentum trans- 
fer rate coefficient values. It can be seen that ignoring the 
thermal contribution to the momentum transfer rate coef- 
ficient decreases it typically by a fac tor of > 6. We also plot 
in this figure the value obtained by Geiss fc Buergi (1986 ) 
illustrating the uncertainties in the momentum transfer 
rate coefficient (more on this in Appendix A. 2). 



5. Concluding remarks 

We performed detailed non-equilibrium calculations of the 
thermal and ionization structure of atomic, self-similar 
magnetically-driven jets from keplerian accretion disks. 
Current dissipation in ion-neutral collisions - ambipo- 
lar diffusion-, was assumed as the major heating source. 
Improvements over the original work of (Safier 1993a 1)) 
include: a) detailed dynamical models by Ferreira (1997) 
where the disk is self-consistently taken into account but 
each magnetic surface assumed isothermal; b) ionization 
evolution for all relevant "heavy atoms" using Mappings 
Ic code; c) radiation cooling by hydrogen lines, recom- 
bination and photoionization heating using Mappings Ic 
code; d) H-H + momentum exchange rates including ther- 
mal contributions; and e ) more detailed dust description. 

|5afier , warm jets with a hot tempera- 
'■ 4 K. Such a plateau is a robust 



We obtain, as 
ture plateau at T ~ 10 4 
property of the atomic disk winds considered here for ac- 
cretion rates less than a few times 1O _5 M0 yr _1 . It is 
a direct consequence of the characteristic behavior of the 



wind function G(x) defined in Eq. gq: (i) G(x) increases 
first and becomes larger than a certain value fixed by 
the minimum ionization fraction (see Fig. |^); (ii) G(x) 
is flat whenever ionization freezing occurs (collimated jet 
region). More generally, we formulate three analytical cri- 
teria that must be met by any MHD wind dominated by 
ambipolar diffusion heating and adiabatic cooling in order 
to converge to a hot temperature plateau. 

The scalings of ionization fractions and temperatures 
in the plateau with M acc and wq found by Safier are re- 
covered. However the ionizations fractions are 10 to 100 
times smaller, due to larger H-H + momentum exchange 
rates (which include the do minant thermal velocity con- 
tribution ignored by Safiei) and to different MHD wind 
dynamics. 

We performed detailed consistency checks for our solu- 
tions and found that local charge neutrality, gas thermal- 
ization, single fluid description and ideal MHD approxi- 
mation are always verified by our solutions. However at 
low accretion rates, for the base of outer wind regions 
(ruo ~ 1AU) and increasingly for higher £, single fluid cal- 
culations become questionable. For the kind of jets under 
study, a multi-component description is necessary for field 
lines anchored after zuo > 1 AU. So far, all jet calculations 
assumed either isothermal or adiabatic magnetic surfaces. 
But our thermal computations showed such an increase 
in jet temperature that thermal pressure gradients might 
become relevant in jet dynamics. We therefore checked 
the "cold" fluid approximation by computing the ratio of 
the thermal pressure gradient to the Lorentz force, along 
(/3||) and perpendicular (/3±) to a magnetic surface. Both 
ratios increase for lower accretion rates and outer wind 
regions. We found that for some solutions, thermal pres- 
sure gradients play indeed a role, however only at the wind 
base (possible acceleration) and in the recollimation zone 
(possible support against recollim ation) . Fo rtunately, (as 



will be seen in a companion paper, Paper II ), the dynam- 
ical solutions which are found inconsistent are also those 
rejected on an observational ground. Therefore, it turns 
out that the models that best fit observations are indeed 
consistent. 
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Appendix A: Multicomponent MHD equations 

A.l. Single fluid description 

Let us consider a fluid composed of a species of numeri- 
cal density n a , mass m Q , charge q a and velocity v a . All 
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species are assumed to be coupled enough so that they 
have the same temperature T. To get a single fluid de- 
scription, we then define 



P 
pv 
P 
J 



Y, a n Q kBT 



(A.l) 



as being the flow density p, velocity v, pressure P and 
current density J. We consider now a fluid composed of 
three species, namely electrons (e), ions (i) and neutrals 
(n) . The equations of motion for each specie are 



Dv, 



Dt 
Dvi 
1 Dt 



Dvc 

: Dt 



r = -VP n - P „V$ 



G 



-VP l - Pi V$ G + F e 



qirn(E H x B) 

c 



VP C - 

-e n c (E 



V. 



- x B) 



(A.2) 



(A.3) 



(A.4) 



where $g is the gravitational potential and the colli- 
sional force of particles a on particles (3 is given by 
F a/3 = m a pn a v a p(v a - vp), m a p = m a mp/{m a + mp) 
being the reduced mass, v a p — np (a a pv) the collisional 
frequency and {(J a pv) the averaged momentum transfer 
rate coefficient. 

A single fluid dynamical description of several species 
is relevant whenever they are efficiently collisionally cou- 
pled, namely if they fulfill ||« a =e,n,t — v \\ "C \\v\\. Under 
this assumption and using Newton's principle (E a ^pF a p = 
0), we get the usual MHD momentum conservation equa- 
tion for one fluid 

p— = -VP-pV<S> G + -J x B (A.5) 
Dt c 

by adding all equations for each specie. The Lorentz force 
acting on the mean flow is 

-Jx B = (l + X)(F in + F en ) , (A.6) 
c 

where X = pij p n . Even if the bulk of the flow is neutral, 
collisions with charged particles give rise to magnetic ef- 
fects. In turn, the magnetic field is coupled to the flow by 
the currents generated there. This feedback is provided by 
the induction equation, which requires the knowledge of 
the local electric field E, Its expression is obtained from 
the electrons momentum equation 

v 1 

E + — x B = (m ie niV ie v ie + m ne n n v ne v ne ) 

c en c 



en c 



(A.7) 



where v a p = v a — vp is the drift velocity between the two 
species. Due to their negligible contribution to the mass 



of the bulk flow, all terms involving the electrons iner- 
tia have been neglected (electrons quite instantly adjust 
themselves to the other forces). 

All drift velocities can be easily obtained. The electron- 
ion drift velocity is directly provided by Vi e — J/en c . 
Using Eq. (A.6) and noting that p c v en <C p%Vi n w e get the 
ion-neutral drift velocity 



-J x B 



(i + x) 



en,-. 



(A.6 



On the same line of thought, the electrons velocity is v e 
v — (v — v Q ) where 



v r . 



X 



:(Vi 



l + X 

J 

en c (1 + X) 2 m m riiV, 



Vr, 



l + X' 

1 -J x B 



(A.9) 
(A.10) 



Gathering these expressions for all drift velocities, we ob- 
tain the generalized Ohm's law 



E 



v „ r ±JxB VP C 
— x B — rjj + £ 

c en e en e 



1 ^(JxB)xB 
(l + X)* 



(A.ll) 



where r\ = (m ne n n v ne + mi e riiVi e ) / (en e ) 2 is the electri- 
cal resistivity due to collisions. The corresponding MHD 
heating rate writes 



MHD 



= J • E' = J ■ E 



— x B 



(A.12) 



where E' is the electrical field in the comoving frame. This 
expression leads to equation (p"7|). 

The generalization of this derivation for a mixture 
of several chemical elements has been done in a quite 
straightforward way. The bulk flow density becomes p = 
7k + J>n + pet where the over line stands for a sum over 
all elements (ions and neutrals), with X = jH/jhl- The 
neutrals and ions velocities are means over all elements, 
(v n ,i) = J2 n ,i Pn,iV n ^jjT~. The conductivity and colli- 
sion terms are also sums over all elements, namely fj = 
(m ne nnv ne +mi e riiv le ) I (en c ) 2 and m ln mvin, and are com- 
puted using the expressions for the collision frequencies. 

A.2. Momentum transfer rate coefficient 

For ion-electron collisions we use the canonical from 
[Bchunk (1975 ), summed over all species: 



27re 4 



3(fc B Tc) 2 



8fc B T 



J2riiZflnAi (A.13) 



with the Coulomb factor A, = (3/2Z,e 3 )^ (k B T c ) 3 /un c . 
For the collisions between electrons and neutrals we 



use the expression of pstcrbrock (196l| ) for the colli- 
sional momentum transfer rate coefficient between a neu- 
tral and a charged particle, which corrects the classical one 
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(eg., Schunk 1975) for strong repulsive forces at close dis- 
tances. Its expression is (av) n ,i—e = 2.4l7re-y/a n /m ni j_ e , 
where the polarizabilities a n used are also taken from 
Ostcrbrock. We thus obtain 



m en n n v ne = 2.417T e n e ^ n ny /m c a r , 



(A.14) 



=H,Hc 



Finally, it is mainly the ion-neutral collision momen- 
tum transfer rate coefficient determines the ambipolar dif- 
fusion heating. It can be computed with the previous mo- 
mentum transfer rate coefficient expression. However as 



noted by [Draine (1980| ) the previous expression under- 
estimates a at high velocities. Thus, as Draine, we take 
the "hard sphere" value for the cross-section (as — 10 -15 
cm 2 ) whenever it is superior to the polarizability one. 
For intermediate to hight ionizations (/ H + > 10 -4 ) the 
dominant ion- neutral collisions are those between H-H + . 
Charge exchange effects between these two species will 
amplify (tr H H + v) above the values expected by polariz- 



ability alone and thus it is computed separately (Eqs. A. 16 



and A. 17). We thus obtain for ion-neutral collisions 



+ 22 niirim max (as v; (<xw)h,i) 



i>H+ 



(A.15) 



+ ^ nimm e max (a s v; (av) He ,<) 



where v — ^J^k^T j Trm,i n + vf n . For v < 1000 km s 1 . 

The value of (chh+ u ) which we used is given by 
Draine (1980| ), 



(°" H H+^) 

1 cm 3 s _1 



3.26 x 10~ 9 
2.0 x 10- 9 v< 



v < 2 km s 
v > 2 km s~ 



(A.16) 



Saficr (1993a) used the expression v = Vi n which, as dis- 
cussed in section 4.8, results in a smaller momentum trans- 



fer rate coefficient. Geiss & Buergi (1986) computed an- 
other expression of the H-H + momentum transfer rate co- 
efficient, which provides 



(a Hn +v) = 1.12 x 10" 8 T 4 2 (l-0.121ogT 4 ) 2 
+ (2.4 -0.34(1 + 2 logT 4 ) 2 ) x 10~ 9 cmV 1 



(A.17) 



In Figure ^ we compared both momentum transfer rate co- 
efficients, they typically differ in 40%, which can be used 
as an estimate of their accuracy. It is thus the uncertainty 
in the H-H + momentum transfer rate coefficient that dom- 
inates the final intrinsic uncertainty of our calculations. 

Appendix B: Dust implementation 



As shown by |Saficr| if there is dust in the disk, the wind 
is powerful enough to drag it along. Thus disk winds are 
dusty winds. Dust is important for the wind thermal struc- 
ture mainly as an opacity source affecting the photoion- 
isation heating at the wind base. To compute the dust 




Tsub=1 400K - silicate 
Tsub=1500K - MRN mix 
Tsub=1 750K - graphite 



Fig. B.l. Dust sublimation surfaces geometry for the 
adopted radiation field. 

opacity we need a description of its size distribution, its 
wavelength dependent absorption cross-section and the in- 
ner dust sublimation surface. In the inner flow zones and 
for high accretion rates the strong stellar and boundary 
layer flux will sublimate the dust, creating a dust free in- 
ner cavity (see figure |b|). Results on the ev olution of dust 
in accretion disks by [Bchmitt et al. (1997 ) show that at 
the disk surface the initial dust distribution isn't much 
affected by coagulation and sedimentation effects. Thus 
we assume a MRN dust distribution (Mathis et al. 1977; 
|Draine fc Lee 1984| ): 



drii = riH Ai 



-3.5 



da 



(B.l) 



where drii is the number of particles of type i ( "astronom- 
ical silicate" - Sil or graphite - C) with sizes in [a, a + da], 
and 0.005/Ltm < a < 0.25/mi, A S u = 10" 2511 cm 2 ' 5 H" 1 
and A c = 10~ 25 - 16 cm 2 5 H _1 . We then proceed by av- 
eraging all relevant grain quantities function of size and 
species (Fi(a)) by the size/species distribution, 

(Fi(a)) a = ]T ^(a)^ (B.2) 

i=Sil,G 



In order to compute the sublimation radius some de- 
scription of the dust temperature must be made. For sim- 
plicity, we assume the dust to be in thermodynamic equi- 
librium with the radiation field, the dominant dust heat- 
ing mechanism. In our case, the central source radiation 
field will dominate throughout the jet, except probably 
in the recollimation zone, where the strong gas emission 
overcomes the central diluted field. However in this region 
dust is no longer relevant for the gas thermodynamics and 
we will therefore only consider dust heating by the central 
source. The dust temperature T gr for a grain of size a is 
obtained by equating the absorbed to the emitted radia- 
tion (eg., |Ticlcns fc Hollcnbach 1985| ), 



Ana 2 (Q a ) (T gr )aT 4 



Qf s (v)AnJ v dv (B.3) 



where a is t he grain size (Q a ) (T gr ) is the Planck-averaged 
emissivity ( Draine & Lee 1984); |Laor & Draine 1993 ; 



Draine & Malhotra 1993), a is the Stefan-Boltzmann con- 
stant, Qf s (v) is the dust absorption efficiency]^] and 4-7T J„ 

2 The dust absorption efficiency is related to the dust ab- 
sorption cross-section by = ■Ka 2 Q'^ aa (v). 
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is the central source radiation flux at the grain position 
given by equation Averaging out the previous equation 
by the size/species distribution (eq. B.2) we obtain, 



4(gr(r gr )X : 



(B.4) 



where we describe the central source flux by F{v) which 
is attenuated only by the dust opacity r. For simplicity 
F{y) is taken as exactly the same as in (Saner , ie. a classi- 
cal boundary layer ( Bcrtout ct al. 1988| ) . The sublimation 
radius is obtained from the previous expression by noting 
that at its position r v = 0, 



sub(fl)) 



l9*(0) (QZ hs (T*))Tt + g hl (8) (Q*»(T bl ))T. 



4(Qr(r sub ))T s 4 ub 



(B.5) 



where Tbi and T* are the boundary layer and star tem- 
peratures, i?* the stellar radius and gh\{8) 1 9*{9) are the 9 
dependent terms of the radiation field (given in Bertout et 
al. and Safier). We assume a dust sublimation temperature 
T sub of 1500 K. 

With the dust sublimation radius in hand we can now 
proceed to compute the dust optical depth defined as, 



Tv{r,6) 



AO) 



Kl bs {r')dr' 



AO) 



nn(r',9)a(v) a dr' 
(B.6) 



where r m (9) is the radius inside which there is no dust. 
This radius is given by the inner flow line r ro . (9) and by 
the sublimation radius (r su b(9)) (see figure such that 
Tin(9) = max( (r su b) ; The dust absorption cross- 

section (a(v) a ) is, 



Tra 2 [Q^ s (a,v)A S n 
+ Q^ s (a,iy)A c ]a~ 3 - 5 da 



(B.7) 



Usin g th e self-similarity oin^r, 9) we can integrate equa- 
tion 



B.6 to obtain, 



7in(0)) 



1 



(B.8) 



which was used in equation [T^. We note that at large 
distances from the source, the optical depth converges to a 
finite value, proportional to M acc and whose 9 variation is 
function of the self similar wind solution and central source 
radiation field. Thus for high accretion rates, although the 
central source radiation hardens, the outer zones of the 
wind base are less photoionized than for smaller accretion 
rates. 

Appendix C: Consistency checks 

C.l. Dynamical assumptions 

First, local charge neutrality is always achieved. For exam- 
ple, we achieve a maximum Debye length of r^, ~ 10 5 cm 



\ [V; V c |/V 



"■ I ' I «T ■ ""m i 
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Fig. C.l. Top left: We plot the relevant drift speeds 
normalized to the fluid poloidal velocity. The worst case 
of one fluid approximation violation is obtained for model 
C, Macc = lO^Moyr- 1 and zu Q = 1 AU. Top right: We 
plot the t Q( 3 versus dynamical Td yn time-scales (s) versus 
X, normalized to the Keplerian period at the line footpoint 
(for a 1 M Q star). We plot only the longer time-scales 
(r in = T n i/fi and T cn = r nc // c ). The worst case is obtained 
again for model C, M acc = 10 _8 M Q yr~ 1 and zjq = 1 AU. 
Bottom left: Ideal MHD tests for the worst situation 
(model C, M acc = 10" 8 M o yr- 1 and w = 1 AU). As 
expected from our heated winds, the ambipolar diffusion 
term is the dominant one. Bottom right: ratios of the 
thermal pressure gradient to the Lorentz force versus x> 
along (/3||, solid) and across dash) a magnetic surface 
anchored at tuq. The worst case for /3ii is obtained for 
model A, zu = 1 AU, M acc = lO _8 M yr _1 , the best for 
model A, m = 0.1 AU and M acc = lO -5 M yr -1 . With 
respect to /3j_, the worst case is for model C, zuq = 1 
AU and M acc = 10 _8 MQyr~ 1 , while the best is for model 
A, txj = 0.1 AU and Af acc = 10~ 5 M©yr _1 . Although 
definitely not negligible in some models, those compatible 
with observations do not show important deviations from 
the "cold" jet approximation. 



at the outer radius of the recollimation zone (model C, 
lowest M acc ). 

Second, single fluid approximation requires that rela- 
tive velocity drifts of all species (a = ions, electrons, neu- 
trals) || u— v a ll/H v\\ are smaller than unity. These drifts are 
higher for lower accretion rates and at the outer wind base 
(due t o th e decrease in density and velocity, see Eq. ^ . In 
figure CA we present the worst case for the drift veloci- 
ties, showing that our jets can be indeed approximated by 
single fluid calculations. 

We assumed gas thermalization, which is achieved only 
if collisional time-scales between species r a p — \/v a p 
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are much smaller than the dynamical time-scale Tdyn = 
'cuo(dv z /dx)~ 1 ■ In the collision network considered here, 
the longer time-scales involve collisions with neutrals. 
However, even in the worst situation (see figure C.l), after 
the wind base they remain comfortably below the above 
dynamical time-scale. 

Our dynamical jet solutions were derived within the 
ideal MHD framework. This assumption requires that all 
terms in the right hand side of the generalized Ohm's law 
(equation A.ll) are negligible when compared to the elec- 
tromotive field v x B/c. We consider Ohm's term ||?yJ||, 
Hall's effect || J X B\\/cen c and the ambipolar diffusion 
term (^-) 2 ||(J x B) x B\\/c 2 m in niV in (effects due to the 
electronic pressure gradient are small compared to the 
Lorentz force — Hall's term — ). In figure C.l we present 
the worst case for our ideal MHD checks. We find that 
deviations from ideal MHD remain negligible, despite the 
presence of ambipolar diffusion. As expected, this is the 
dominant diffusion process in our (non turbulent) MHD 
jets. Ambipolar diffusion is larger for low accretion rates 
and at the outer wind base (because the ratio of the am- 
bipolar to the electromotive term scales as (M acc /j) -1 . 

The worst case for the previous three tests is, 
as expected, for the model that attains the low- 
est density: Model C, with the lowest accretion rate 
{Macc = 10~ 8 MQyr^ 1 ) and at the outer edge footpoint 
(nj = 1 AU). 

The dynamical jet evolution was calculated under the 
additional assumption of negligible thermal pressure gra- 
dient (cold jets). Since it is the gradient that provides a 
force, one should not just measure (along one field line) 
the relative importance of the gas pressure to the mag- 
netic pressure (usual (3 = P/ (B 2 /8ir) parameter). Instead, 
we compare the thermal pressure gradient to the Lorentz 
force, along (/3||) and perpendicular (f3±) to the flow, 
namely 



01 



F \\ 

F± 



VP 



v p ■ (J x B) 

Va • VP 
'Va • (J x B) 



(CI) 
(C.2) 



Here a(tu, z) is the poloidal magnetic flux function, hence 
Va is perpendicular to a magnetic surface. High values of 
imply that the thermal pressure gradient plays a role 
in gas acceleration, whereas high values of f3± show that 
it affects the gas collimation. 



In figure C.l we plot the worst case of cold fluid vio- 
lation and best case of cold fluid validity. Again the worst 
case appears at lower accretion rates and in outer wind 
zones. It can be seen that high values of (3± and 0\\ can 
be attained, hinting at the importance of gas heating on 
jet dynamics (providing both enthalpy at the base of the 
jet and/or pressure support against recollimation further 
out). We underline that models inconsistent with the cold 
fluid approximation are those found to have t he largest dif- 
ficulty in meeting the observations ( Paper II ). Conversely, 
models that better reproduce observations also fulfill the 
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Fig. C.2. Ignored heating/cooling terms (in erg 
s _1 cm -3 ). Left: We compare the cooling term 

3 l,~r 



~khTDf c /Dt (dashed) with adiabatic cooling A a d 



la 



(solid) for the worst case (model C, M acc = lO _8 M0yr _1 
and ruo = 0.1 AU). Right: Grain heating/cooling 
rate |r gr | (dashed) compared with ambipolar diffusion 
heating Tdrag (solid) for the worst case (model A, 
M acc = lCT'^Moyr- 1 and vo = 1 AU). T gr is only 
significant at the very base of the wind, and will not 
affect the thermal state further out in the jet. 

cold fluid approximation. For those models, the thermal 
pressure gradient appears to be fairly negligible with re- 
spect to the Lorentz force. 

C.2. Thermal assumptions 

Finally, we check that all ignored heating/cooling pro- 
cesses are not relevant when compared to adiabatic cooling 
and ambipolar diffusion heating. 

The first ignored process is the term —^khTDf e /Dt. 
This term decreases for increasing accretion rate and zuq 
due to the lower ionizations found in these regions. It 
is plotted in figure |C.2| for the worst case (model C, 
Macc = 10 _8 M Q yr -1 and m = 0.1 AU). There, it reaches 
at most ~13% of A a di a . Typical values for higher accretion 
rates are only ~ 0.1% of A a di a . 

Next we consider heating/cooling of the gas by col- 
lision with dust grains, given by Hollenbach & McKee 
1(19791 ): 



r gr = n {n„a„) 



8k B T 

7T771H 



/(2fc B T g r - 2k B T) (C.3) 



where (n gI a gI ) is computed from the adopted MRN dis- 
tribution, and / = 0.16 is the sticking parameter that 
takes into account charge and accommodation effects for 



a warm gas ( Hollenbach fc McKee 1979 ). With these val- 



ues the grain heating/cooling becomes 
r gr = 4.78 x 10 _34 n 2 (T gr - T) 



erg s 



(C.4) 



This ter m inc reases with increasing vuq and accretion rate. 
In figure |C.2] we compare it (in absolute value) with Tdrag 
in the case where its contribution is the most important 
(model A, M acc = lO _5 M yr _1 and zu Q = 1 AU) 



r gr is 
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initially positive (dust hotter than the gas), but changes 
sign at x — 0-4, where the gas becomes hotter than the 
dust, becoming an effective cooling term. It is only sig- 
nificant at the very base of the wind, where it exceeds 
the ambipolar diffusion heating term by a factor ~ 3.5. 
However, this effect will not have important consequences 
in terms of observational predictions: We will show in next 
section that the thermal state in the hot plateau (where 
forbidden line emission is excited) is not sensitive to the 
initial temperature. Furthermore, the outer streamlines at 
zuq ~ 1 AU contribute much less to the line emission than 
the inner ones. At lower accretion rates < 10 _6 M©yr _1 , 
r gr is always <10 % of Tdrag- 

Heating due to cosmic rays, which could be impor- 
tant i n the out er tenuous zones of the wind is ( [Spitzer "fe 
Tomas ko 196^) , 

T cr = n H (AE = 1.9 x Kn 28 n# erg s" 1 cm~ 3 (C.5) 

where £ is the ionization rate which we took as £ = 3.5 x 
lCT 17 s" 1 ( |Webber 199§| ) and AE = 3.4 eV is the average 
t hermal energy transmitte d to the gas by each ionization 
flSpitzcr fe Tomasko 19681) . This effect is at most ~ 3.6 x 
10~ 6 times r dr ag for model A, M acc = 10 _5 M Q yr _1 and 
vj = 0.1 AU. 

Finally the thermal conductivity along magnetic field 
lines was computed with the Spitzer conductivity for a 
fully ionized gas ( Lang 1999 ) and is irrelevant (at most 
~ 10~ 6 of the adiabatic cooling term) the maximum be- 
ing achieved at the recollimation zone where the physical 
validity of our MHD solutions ends. It should be pointed 
out that Nowak fc Ulmschneider (1977 ) compute the ther- 
mal conductivity for a partially ionized mixture in ion- 
ization equilibrium and found that for low temperatures 
T ~ 10 3_4 K the Spitzer expression underestimates the 
conductivity by a factor 10 2 . However this is still too small 
to be important. 

C.3. Dependence on initial conditions 

Formally, our temperature integration is an initial value 
problem. In the absence of a self-consistent description 
of the disc thermodynamics, there is some freedom in the 
initial temperature determination. It is therefore crucial to 
check that the subsequent thermal evolution of the wind 
does not depend critically on the adopted initial value. 

Satier obtained the initial temperature by assuming 
the poloidal velocity at the slow magnetosonic point 
(i> PjS ) to be the sound speed for adiabatic perturbations 
T s = /imnVp s /jk-B- Here, we have chosen to compute the 
initial temperature assuming local thermal equilibrium 
DT/Dt = 0. Our method produces lower initial temper- 



©„ = 0.1 AU 



O n = 1 AU 



atures than Safier due to adiabatic cooling. 

For high accretion rates M acc > lO _6 M0yr _1 our ini- 
tial temperature versus wq has a minimum at the begin- 
ning of the dusty zone: Inside the sublimation cavity, the 
thermal equilibrium is between photoionization heating 
and adiabatic cooling. Just beyond the dust sublimation 
radius, photoionization heating is strongly reduced, but 
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Fig. C.3. Effect of initial temperature on the thermal evo- 
lution for model B, Af acc = lO _6 M0yr _1 . The several 
initial temperatures are in dashed 50 K, 100 K, 500 K, 
1000 K, 2000 K and 3000 K. In solid we plot the solution 
obtained by our standard initial conditions. Note that the 
almost vertical evolution of the temperature for very low 
initial temperatures is not an artifact. The field line an- 
chored at zuq = 0.1 AU crosses the sublimation surface at 
X ~0.5. 



the ionization fraction is still too high for efficient drag 
heating, resulting in a low initial equilibrium temperature. 

The initial ionization fraction is similarly determined 
by assuming local ionization equilibrium Df\/Dt = 
for all elements. It decreases with zuo. For M acc = 
lO _5 M0yr _1 and wq > 0.8 AU, the initial ionization frac- 
tion is set to a minimum value by assuming that Na is fully 
ionized, which is somewhat arbitrary. However, as gas is 
lifted up above the disk plane, the dust opacity decreases 
and the gas heats up, so that ionization becomes domi- 
nated by other photoionized heavy species and by protons, 
all computed self consistently. 

In order to check that our results do not depend on 
the initial temperature, we have run model B for a broad 
range in initial temperatures. As shown in Figure |C.3[ ), 
we find that the thermal and ionization evolution quickly 
becomes insensitive to the initial temperature. If we start 
with a temperature lower than the local isothermal condi- 
tion, the dominant adiabatic cooling is strongly reduced, 
and the gas strongly heats up, quickly converging to our 
nominal curve. If we start with a higher initial tempera- 
ture, adiabatic cooling is stronger, and we have the charac- 
teristic dip in the temperature found by Safier . Our choice 
of initial temperature has the advantage o f red ucing this 
dip, which is somewhat artificial (see figure 03). In either 
case, we conclude that our results are robust with respect 
to the choice of initial temperature. In particular, the dis- 
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tance at which the hot plateau is reached, which has a 
crucial effect on line profile predictions, is unaffected. 
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